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The interplay of quantum statistics, interactions and temperature is studied within the framework 
of the bosonic two-component theory with repulsive delta-function interaction in one dimension. We 
numerically solve the thermodynamic Bethe Ansatz and obtain the equation of state as a function 
of temperature and of the interaction strength, the relative chemical potential and either the total 
chemical potential or a fixed number of particles, allowing to quantify the full crossover behaviour 
of the system between its low-temperature ferromagnetic and high-temperature unpolarized regime, 
and from the low coupling decoherent regime to the fermionization regime at high interaction. 



I. INTRODUCTION 

The increasingly common experimental realization of 
interacting quantum systems using cold atoms has re- 
cently reignited interest in pushing our understanding 
of many-body physics beyond the traditional mean-field 
level [T] . This aspect is most in prominence in effectively 
one-dimensional realizations of bosonic ^^Rb quantum 
gases with tunable interaction strength |lHS]j for which 
the whole crossover from weakly- to strongly-interacting 
physics is accessible. 

The case of locally interacting atoms confined to a uni- 
form one-dimensional channel has the theoretical pecu- 
liarity of being integrable. The simplest case of a sin- 
gle bosonic species defines the well-known Lieb-Liniger 
model [71 IH] for which a recent experimental study has 
shown that the observed thermodynamic properties can 
be understood from the theory of integrable systems at 
finite temperatures |9j. In order to study even richer 
highly correlated systems, multicomponent (spinor) gases 
have been experimentally realized [TUHH]. This exten- 
sion involves different hyperfine states which provide a 
pseudospin degree of freedom [T3HI6] . The control of the 
intra- and inter-species interaction strength via Feshbach 
resonances or state-dependent potentials P^THTO] opens 
the way for realizing a variety of integrable models. 

The main interest in pursuing the study of multicom- 
ponent systems is that they provide situations where im- 
portant interaction and quantum statistics effects coex- 
ist. These two aspects are not unrelated even in single- 
component systems: as a simple illustration, for a sin- 
gle bosonic species in one dimension, the limit of in- 
finitely strong interactions (impenetrable bosons) causes 
a crossover from bosonic to effectively fermionic be- 
haviour [201 m] , at least for physical quantities of density 
type. Considering more than one component however 
opens the door to much richer effects like the presence of 
spin wave excitations, with the possibility of crossover to 
many more regimes than the one component case. 

The study of multicomponent integrable systems re- 
ally begins with the spin- 1/2 fermion problem p2Vl26j . 
An interesting feature of this system is that the attrac- 



tive case has a correspondence to a bosonic system with 
twice the interaction, in the sense that the equation for 
the ground state and particle energy coincide up to a 
sign |26j . A fundamental step forward was thereafter 
achieved by Yang, who showed that the repulsive delta- 
function interactioir problem admitted an exact solution 
irrespective of the symmetry requirement imposed on the 
wavefunction [57]. For spin- 1/2 particles, he showed that 
a generalized Bethe hypothesis in the form of what is to- 
day called a nested Bethe ansatz could provide the sys- 
tem's wavefunctions, obtained the continuum equations 
for the ground state, and calculated the general bound- 
state S'-matrices [55] ■ Sutherland [29] generalized this to 
any irreducible representation of the permutation group, 
so that in particular systems of type (in his notation) 
B^py with X species of bosons and y species of fermions 
were amenable to an exact solution [30 . The ground 
state and excitations of multicomponent fermionic sys- 
tem were studied, both for repulsive and attractive in- 
teractions, by Schlottmann [311 131] • This made exten- 
sive use of the 'string hypothesis' for the solutions to the 
Bethe equations, yielding the various dispersion branches 
in the repulsive case, and the gapped color singlet ground 
states in the attractive one. 

The bosonic multicomponent case has been less exten- 
sively studied. For two components and in contrast with 
the Fermi gas, the ground state is (pseudospin) polarized 
as expected from basic arguments f30j or more formally 
from a general theorem valid when spin-dependent forces 
are absent [33]. In pseudospin language, the ground state 
is thus ferromagnetic, and the excitations at large cou- 
pling correspond to those in an isotropic XXX ferro- 
magnetic chain j34j , revealing thermodynamic properties 
which are drastically different from those of the one- 
component Lieb-Liniger gas [35] |3B]. Furthermore, re- 
cent studies of the ferromagnetic ground state revealed 
different dispersions for the charge and pseudospin exci- 
tations which therefore exhibits a spin-charge separation 
[371 [35] . The regime of strong interaction is still not 
completely understood and even if Girardeau's Fermi- 
Bose mapping has been showed and used for the study of 
the ID spinor Bose gases [33] leading to a paramagnetic 
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Tonk-Girardeau regime [31], this approach fails to pro- 
vide first correction terms to this fermionization regime 
since the discernabihty of the bosons is missing. 

The purpose of our paper is to further stimulate the 
contact between integrability theory and experiments on 
multicomponent cold atoms, by providing quantitative 
predictions for the equation of state and population den- 
sities of two-component interacting Bose gases as a func- 
tion of temperature, interaction strength and of either 
the available chemical potentials or the chemical poten- 
tial difference and a fixed density of particles. This work 
broadens and extends our earlier paper |40j . The pa- 
per is organized as follows. In section [llj after defining 
our notations, we quickly review the construction of the 
eigenstates of the theory using the Bethe Ansatz, and 
how these can be used to obtain the thermodynamics 
of the system in the continuum limit via the solution 
of a (infinite) set of coupled integral equations. Section 
|III| outlines the method we have used to solve this sys- 
tem numerically, using two different approaches allow- 
ing cross-checking of the results. Section IV discusses 



the effect of the thermal fluctuation over the ferromag- 
netic ground-state, whereas section [V| provides results 
on the more challenging intermediate regimes. Section 
VI discusses the results in the decoherents regime of low 
coupling where we compare the numerical results with a 
perturbative result and section [VlI| prescnts the results at 
strong coupling where the gas enters the fermionization 
regime. We end with conclusions and perspectives. 



II. SETUP 

Consider a collection of bosonic atoms of equal mass 
but having an internal SU (2) degree of freedom (in prac- 
tice, this would be for example two distinguishable hy- 
perflne states, and could be thought of as a (pseudo- 
)spin— 1/2. The unique feature differentiating from a sin- 
gle species is the fact that distinguishability imposes sym- 
metry of the many-particle wavefunction on the same- 
spin particles only. For definiteness, we consider a one- 
dimensional ring of length L, in which a total of N atoms 
circulate. The first-quantized Hamiltonian of the system 
includes a free dynamical term to which a spin-blind in- 
teraction term is added, and reads 
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The effective one-dimensional coupling parameter gio is 
related to the effective ID scattering length aiu [ill via 
the relation gm = lT?aio /^m. Hereafter, we will use the 
effective interaction parameter c = giDrn/h^, and adopt 
the traditional convention of setting h = 2m = 1 to sim- 
plify the notations. Note that this choice of interaction 
term involves fine-tuning two parameters: more gener- 
ally, we could have different intra- and inter-species scat- 
tering lengths. To preserve integrability however, these 



must all be equal. 

Specializing to N atoms of which Af have (in the 
adopted cataloguing) spin down, the Bethe Ansatz pro- 
vides eigenfunctions fully characterized by sets of rapidi- 
ties (quasi-momenta) kj, j = 1,...,N and pseudospin ra- 
pidities Xa, a = 1, M, provided these obey the N + M 
coupled equations p71 [53] 
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for j = 1,...,N and a = 1,...,M. For a generic eigen- 
state, the solution to the Bethe equations is rather in- 
volved and cannot be obtained in closed form. Two ob- 
servations allow to push the treatment further: first, for 
c > 0, the kj rapidities live on the real axis. This is not 
true of the Aq which are found to be generically complex, 
but arranged into regular patterns called strings [131 [33] ■ 
An n-string of A's is a congregation of n rapidities sharing 
the same real value and having an even spacing of height 
c in the imaginary direction. The adopted notation for 
the a-th member of a n-string labeled by the index a and 
centered on A^ is thus A^'° = AJJ -|- if (n -|- 1 - 2a), with 
the equality being exact only up to deviations which (ac- 
cording to the traditional string hypothesis) vanish in the 
infinite size limit. Throughout our work, we will adopt 
this as a working hypothesis. Since the total number of 
each type of string is conserved under time evolution, 
each string type represents a quasiparticle of the theory. 
In the thermodynamic limit, the distribution of all ra- 
pidities can be encoded into a set of smooth functions 
representing the densities of roots for each string type. 
The Bethe equations then become a set of coupled in- 
tegral equations for (quasi)particle and (quasi)hole root 
distribution functions. We refer the reader who is unfa- 
miliar with these to our summary of important formulas 
in appendix ( [A] ). 

The Thermodynamic Bethe Ansatz (TBA) allows to 
exploit the condition of thermal equilibrium O [JS] to 
obtain the Yang- Yang- Takahashi (YYT) like equations 
k34ji35ji45j for e(A), the dressed energy, and e„(A;), length- 
n string dressed energy, n — 1,2, ... 
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with the standard convolution notation g * h{k) 
J^^dk'g{k — k')h{k'), and the kernels a„(fc) 
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The set of coupled 



equations is completed with the asymptotic conditions 
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for high-level functions. From these can be derived the 
large-rapidity asymptotes 
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for the large rapidity asymptotic values of the individual 
functions, where we have defined the numbers 
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The thermodynamics of the system is provided by the 
solution set of dressed energies as a function of the tem- 
perature T, the total fj,{= ^'-'^^^ with fii the chemical 
potential specific to the ith component) and relative 
(= ^'-"'^^ ) chemical potential (see appendix[c|concerning 
the c parameter). The Gibbs free energy per unit length 
is given by 
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while the linear density of the ith boson component is 
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The entropy density is given by the standard thermody- 
namic identity 
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and the local density-density correlator jJSl HT] is given 
(using the Hellman-Feynman theorem) by 
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III. NUMERICAL TREATMENT 

To solve the infinite system of transcendental coupled 
equations (3]), we have developed two different numer- 
ical algorithms. We can then independently check the 
results by comparison. We first discuss the common ap- 
proach for the numerical treatment. Afterwards, we will 
describe and motivate the choices we made to build the 
two algorithms. 

Two cutoffs are applied on the system ([s]) to implement 
a numerical process. Firstly we reduce the number of 
functions for computing to nmax, replacing for n > rimax 
these functions by their asymptotic value (limit n — > oo 



in (|4|)). Secondly, we reduce the range of integration 
of the convolution to a finite value. Supposing that as 
k -> ±oo, e{k) ^ fc^ and that e„ becomes £^ ([sj, we 
limit the range of the k values to [— A„ : A„]. In con- 
sequence, we estimate the solution by l-l-nmax functions 
{e, ei, £2, . . . , £Tin,ax} evaluate them by Ni points 

over the interval [— : A^] (i = 0, . . . ,ni„ax)- To com- 
pute the solution we proceed by iterations, starting from 
the free bosons form of e and the asymptotic values . 

The previous paragraph describes how to get the par- 
ticle dressed energy, e{k), and consequently allows one 
to compute the Gibbs free energy ([t]). But in order 
to compute any other thermodynamic quantity involv- 
ing a derivative of 5 (|8j |9) [lO|, we use the system for 
the corresponding set of functions | a^^, ir^, . . .1 with 



e {/i,f^,T,c} (Bl 



|B2[ |B3[ ). We achieve the numeri- 
cal solutions by the same method of discretisation intro- 
duced before and using the set {£, £1, . . . , £n^.^^} solution 
of ([3]). It would be possible to numerically differentiate 
G or £ to compute these quantities. This method would 
however achieve only much less accuracy for given com- 
putational effort. 



1 . Fixed density of particles 

For physical interpretation of the results and for the 
identification of the different regime crossovers, results 
with fixed density of particles could be more convenient. 
For this purpose, we implemented a Newton's method 
on top of our main algorithms that finds the chemical 
potential, /i, corresponding to a desired density, ni +n2. 



2. Accuracy and precision 

By the use of the method mentioned above to solve the 
system, we are confronted with two limitations on exacti- 
tude. The first comes form the numerical approximation 
of the system: discretization of the functions, limitation 
of the integration range and the number of function leads 
to an imprecision on the results. Second, the fact that we 
solve the system by iteration approaching but not reach- 
ing the solution, leads to a inaccuracy. 

The imprecision is limited as described hereafter. The 
discretisation of the functions adds an error 0[-^) in 
the convolutions {N being the numbers of points) which 
induces an imprecision that one can easily keep low. Con- 
cerning the range cutoff, all the thermodynamic quanti- 
ties are c>c e"^'*^-'/^, with e{k) ^ fc^ when fc ^ 1. There- 
fore we reduce this effect on the results by taking an 
appropriate k spaci ng such tha t A^/T ^ 1. Finally, as 
we can see from ([3) |B1|B2|B3[ ), all the contributions of 
the £„ functions in £(fc) are oc e"^"^*^-'/"^. Knowing that 
for n 3> 1, Snik) — 2nfl, we therefore keep this impreci- 

— 3> 1. In the results shown 



sion small by choosing — 
in all plots, the precision of the results is estimated to 
always be much smaller than the width of the curves. 
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The problem is different for the accuracy. Supposing 
that the solutions {e,ei,e2, ■ ■ •}, {g^, • ■} (^^^r G 

{fi,fl,T,c}) exist, we suppose that by iteration we ap- 
proach these solutions but the distance from the solution 
is nevertheless unknown. We estimate the accuracy em- 
pirically. Increasing the total number of iterations expo- 
nentially, we observe the convergence of the results and 
judge the accuracy value. In the following results, the ac- 
curacy is estimated to be of order of the line width and 
therefore globally is the higher limitation on exactitude 
of the results. 



A. FFT based algorithm 
1. Idea 

During the iterative process, the major part of cal- 
culation time is taken by the evaluation of each con- 
volution. Indeed by calculating the integrals by simple 
trapezoidal sums, this charge represents ^ 0{NiNi±i) 
operations for each function. Starting from this ob- 
servation, the basic idea of this algorithm is to use 



the Fast Fourier Transform to calculate the convolution: 
(/ * g)(x) = FT-\FT(/) • FT{g)) (with FT: Fourier 
Transform) . The computing time of each convolution us- 
ing FFT is thus only - 0(7V^ \og{N,)). The conditions of 
use of the FFT are that the functions / and g must be 
integrable and that the values of these functions must be 
zero outside [— : A^]. This could be easily achieved by 
treating the constant part of the function separately from 
the nontrivial part. Moreover this method imposes the 
numerical constraints that all the points must be equally 
spaced and that the range and the numbers of points of 
each function must be the same. 



2. Practically 

As a consequence of this, we have a set of Jimax + 1 
functions each to be evaluated on N points and within 
the range [—A, A]. We start by setting up the system 
with three arbitrary parameters: (Dq, Aq, n^j^^) with 
D = i^. During the convergence, we adjust them dy- 
namically with the use of the following precision indica- 
tors. We firstly estimate the nth iteration precision with 
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with e^"'\ki) being the value of e{ki) after n iterations. This formula has to be understood as the variation of the 
Gibbs free energy between two steps (see ([t])). We measure similary how the parameters Umax and A influence the 
precision with these two indicators 
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During computation, if crn„,^^ > cfh, fi-max is increased 
and if ajj > crit,A is lengthened. The precision related 
by the density of points, D, is hard to quantify but can 
be estimated by cross-checking with the second method 
which has a non-uniform distribution of points. We then 
increase it step by step as one goes along the iterating 
process. 
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Since the arbitrary ' with D are reached, we con- 
sider that we have a good evaluation of the solutions. 



Once satisfactory values for af°^,D^°^ are achieved, we 
assume the solution has been reached and we calculate 
the Gibbs free energy from ([t]). The derivatives of the 
Gibbs free ene r gy a re computed using the derivative sys- 



tems (Bl 



B2 



B3) with the final (A, n^ax, -0"°') deter- 



mined by the first iterative process and we approximate 



B. Flexible-density method 

A second, completely independent implementation of 
the numerical solution to the coupled integral equations 
has been pursued as part of our work. Here, we do not 
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make use of the fast Fourier transform, but rather main- 
tain total flexibility in the choice 1) of density of sampling 
points within each function, 2) of the limits used at 
each level, 3) of the relative total number of points used 
at each level, and 4) of the total number of functions 
used. This advantage allows one to concentrate compu- 
tational resources where they are needed, but comes at 
the cost of being only able to perform convolutions be- 
tween e.g. levels i and j at speed of order NiNj where 
Ni is the number of points used at level i. This second 
algorithm performs more or less equally well as the first, 
and allows one to certify the results obtained. 

In summary, this second algorithm works as follows. 
Depending on the physical parameters requested, an ini- 
tial choice is made of the number rtmax of functions to 
be considered, and of the limits at each level. A dy- 
namical parameter called the running precision is initial- 
ized, which estimates the numerical accuracy obtained in 
computing the free energy using the points configuration 
used. The coupled equations are then iterated (possibly 
using extrapolations) in order to achieve a certain degree 
of convergence, measured by the condition that the total 
rate of fiow of all points as the iterations proceed becomes 
smaller than the running precision. 

At this point, a cycle is initiated. This entails a num- 
ber of steps, with the objective of increasing the accu- 
racy, i.e. of decreasing the running precision achieved. 
First, each function is examined in turn, and points are 
added in regions with larger curvature. Second, the lim- 
its Ai are extended (and points added) if the value of the 
function at the previous limit is not sufficiently close to 
its analytically-determined asymptotic value e°°. Third, 
new functions are added (i.e. Umax is increased) if the 
highest function is not sufficiently close to its asymptotic 
value throughout the k line. A new value of the running 
precision is then determined, based on the refinements 
just performed on the distribution of points. Finally, it- 
erations are performed until the fiow rate drops below 
this running precision. 

For a specific set of physical parameters, a total al- 
lowed time is also given to the program. This second 
implementation then performs cycles one after the other, 
yielding increasingly accurate results, until the allowed 
time is exhausted. An estimate of the absolute accuracy 
of the whole procedure can thus be obtained by com- 
paring the results from runs with different total allowed 
times. 



IV. QUANTUM STATISTICS VERSUS 
TEMPERATURE FLUCTUATIONS 

The SU{2) degree of freedom in combination with the 
bosonic statistics in ID leads to a macroscopic behavior: 
the polarization of the ground state. This phenomenon 
which occurs in every bosonic system with no explicit 
component-dependent forces, has been already proven in 
the literature [331 . We will give here an interpretation in 



terms of the string structure of the Bethe solutions and 
provide a quantitative result for the persistence of this 
effect for non-zero temperature. 

The polarization at zero temperature can be directly 
linked to the underlying string structure of the Bethe 
equation solutions. In eq. ([3]), the e(fc) function depicts 
the charge degree of freedom and the £„(fc) functions the 
spin degrees of freedom. Moreover those latter functions 
express the dynamics of quasiparticle forming a colour 2 
energetically disfavored state made of n particles. As the 
temperature of the system goes down, the contribution 
of these states in the equilibrium decreases. 

We hereafter explain how at the limit T = there 
only remains one spin-gapped state gathering all par- 
ticles in the 1st component. In the YYT equations 
([3| where the £„(fc) are the dressed energies of an n- 
string, a phenomenological approach to T = is possi- 
ble. Taking the first line of (|3| and approximating the 
values en{k) « 2nn , if one takes the limit T — >■ 0, 
T'Er^iln [1 +exp(-2^^)] ^ as O is defined posi- 
tive. The contribution of the colour 2 particles then 
disappears from the thermalized state and this reveals 
that the colour 1 component drives away all colour 2 
particles, forming a fully polarized spin-gapped state. 
(In the Bethe equations, only the 2nd component part 
of the wave function is represented by quasiparticles) . 
We show this expelling in figure [l] where we plot polar- 
ization curves (bottom set) as a function of fi. As the 
chemical potential increases, the interaction parameter, 
7 (= „^^„^ ) decreases monotically and the polarization 
persists to higher temperatures. In [38], Fuchs et al. re- 
vealed that the effective mass of an isospin wave above 
the polarized ground state is very high in the strong cou- 
pling regime. Furthermore, it is surprising to see that 
when fi increases, the polarized ground state is more re- 
sistant to thermal fluctuation even though the isospin 
wave mass decreases. 

In the context of spontaneous imbalance in binary mix- 
tures [ISHSOj . it has been shown that at zero tempera- 
ture, a mixed gas is unstable and exhibits a spatial phase 
separation. But no quantitative predictions have been 
made for finite temperature in ID. From figure [T] we see 
that the polarization remains at higher temperatures for 
high value of /i. By qualitative identification of the ferro- 
magnetic behavior with the the spatial demixing, we can 
speculate that a phase separation would resist better to 
temperature in the low coupling regime than for 7 3> 1. 

The figure [T] shows as well curves of polarization for 
fixed particle density, interaction strength and fixed f2 as 
a function of the reduced temperature, t = where 

Td = {ni + is the degeneracy temperature [15] . 
We compare these results with the polarization of an 

interaction- free 2CBG which is ^o_|_„o — t&'ah{D. / k bT) . 
The non-negligeable difference which appears between 
the curves of same $7 is the ferromagntic effect which 
is a consequence of the bosons interaction. 
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FIG. 1: The top graph shows polarization of the 2CBG , 
(ni — n-2)l(ri\ + 712), for fixed density of particles and inter- 
action strength (7) as a function of the reduced temperature 
(r). The set of curves for several different are compared to 
the curves of a free gas. The bottom plot shows the isobar 
polarization as a function of the temperature for four different 
chemical potentials: {—100, 0, 100, 200}. In this latter graph, 
the density 7 and r vary along the curves. 



As T approaches zero, the gas polarizes and the remain- 
ing component behaves like a Lieb-Liniger gas of chemical 
potential /zi = /i -I- fi. The results are then comparable 
to the results of V. N. Popov [5T] for the density of the 
Lieb-Liniger gas at T = and 7 ^ 1: 
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(15) 



As shown in Fig. [2] by lowering down the reduced tem- 
perature of the gas to r <C 1, the total chemical po- 
tential corresponds to the zero temperature low coupling 
regime formula. Moreover, in this particular regime, the 
corresponding polarizations lines turn out to be almost 
constant along 7. 



The specific heat capacity of the gas at fixed density of 
particles which is accessible via the entropy ([9| provides 
a view on the thermal degrees of freedom of the system. 
Figure [3] shows that in a 2CBG at low temperature with 
strong r2, the heat capacity is similar to that of a Lieb- 
Liniger Bose gas contributed by phonons. If the relative 
chemical potential is lower or of order of the tempera- 
ture, the two component degree of freedom appears and 
creates peaks similarly to a paramagnetic spinor Bose gas 
|34j . The maximum in the specific heat moves in higher 
temperature as the relative chemical potential increases. 
The higher temperature results are shown in |4] where we 
can observe that the peaks are located when ^ ^ T . At 
high temperature the gase becomes decoherent classical 



FIG. 2: For 7 < 1, and r < 1, the 2CBG chemical potential 
follows the Popov's expression for the zero temperature Lieb- 
Liniger case. 



(see VI) and the specific heat converges to 1/2 which is 



the value of the simple ID ideal gas. 

Experiments trapping Helium— 4 fluid into ID 
nanopores [52j provide a possible realization of the ID 
Bose gas and give access to measurement of the heat 
capacity of the unidimensional system. A similar realiza- 
tion with an isospin-1/2 might be possible and provide 
measurement of the 2CBG heat capacity. 



y = 1 , n.|+n2 = 1 ,£2 = 0.1 
£1 = 1 
£J = 10 




FIG. 3: The specific heat of the 2CBG for fixed density of 
particles and interaction strength as a function of the reduced 
temperature. The different chemical potentials of each curve 
have a value closed to 7. Similarly to a free 2CBG case, we 
see a peak in the specific heat whos the position depends on 
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FIG. 4: The specific lieat of tiie 2CBG for fixed density of 
particles and interaction strengtli as a function of the reduced 
temperature. The different chemical potentials have a value 
n ^ 7. In the limit of high temperature, the gas becomes 
ideal and the specific heat takes the value 1/2. 



V. RESULTS IN THE INTERMEDIATE 
REGIME 



In this section we will present the results that don't 
belong to a limit regime. Furthermore in those parame- 
ter ranges, we give numerical results for the polarization 
of the 2CBG and the local pair correlator where neither 
the thermal, the charge nor the phase fluctuations dom- 
inate. They compete in the 2CBG state and therefore 
no perturbative approach but only the thermodynamic 
Bethe ansatz can predict a solution. We will discuss and 
describe as much as possible the changes in behaviour 
that occur in these intermediate regimes. 

The set of following graphs in figures [5] and [6] show the 
behaviour of the two component Bose g function of 

7 and any fixed value of i7,T. The values of and T are 
chosen such that in each case the ratio p goes from less 
to more than 1. Following the qualitative description of 
[3S1 ST] for the single-component case, the regimes of the 
gas are identified by the two dimensionless parameters 
7 = — J — and r = 7 — f — ^2 , respectively the interaction 
strength and the reduced temperature. As results pre- 
sented hereafter are made for fixed interaction parameter 
(c), the ratio ^ = ^ is then constant and the regime 
is identified by the position on the 7 axis. At the lowest 
value of 7, 7 < T ^ 1 and the gas quasicondenses in a 
Gross-Pitaevskii (GP) regime with thermal fluctuations. 
At the other end, where 7 > 1, the regime is decoherent 
classical (DC) with r ^ max{l, 7^}. In the case of a qua- 
sicondensate, we see progessively the ferromagntic effect 



with a completely polarized gas whereas the polarization 
reaches the value of an ideal paramagnetic gas when the 
2CBG becomes DC. From this simple view of the data, 
we can try to see how the temperature and the relative 
chemical potential modify these phenomena. 

The first column of figure [5] shows the effects of the 
temperature on the polarization. In the region where 
7 > I the linear density of each component is almost 
classical and the asym ptotic value of the curves are given 
by e;3M!+e'^M2 (^^^ ^I'- Here the charge and coherent fluc- 
tuations are large and hence the statistics and the inter- 
action of the gas don't play any role (the observables de- 
pend only on the temperature and chemical potentials). 
In contrast, for 7 ^ 1 the gas quasicondenses and the 
charge fluctuations vanish. The ratio r/7^ being large, 
the temperature fluctuations exceed the phase fluctua- 
tions and we see that T doesn't influence the polarization 
much. 

The second column of data shows the variations of 
the polarization as a function of fl. The spontaneous 
ferromagnetism in the presence of the quasicondensate 
happens in high interaction strength when the relative 
chemical potential increases. In the YYT equations 
([3]), the effect of Q on the strings appears through the 
asymptotic value of the contribution of the n-strings: 
rX;r=iln [1 +exp(-2^^)]. When fl increases, the 
colour 2 spin-gapped state effect are suppressed and the 
polarization resists higher charge fluctuation (higher 7). 

Figure [6] shows the local density-density correlation 
function as a function of 7 for different values of the rel- 
ative chemical potential and temperature. On the top 
graph the ratio T/7^ is fixed; for 7 ^ 1, the 2CBG is 
thus in a quasicondensate with important thermal fluctu- 
ations. In this regime the gas is ferromagnetic and has 
no effect on the correlation. On the other hand, for large 
values of 7, the gas is DC and the asymptotic value of g'^-* 
follows from Wick's theorem and the Boltzman distribu- 
tion, q\ = 1 + , ' „ — The first order corrections 

" (Eje-^^O 

will be calculated later (VI). In the bottom figure, the 



curves for different temperatures show the nonmonotonic 
behaviour discussed in [iD] . Close to the quasicondensate 
regime, the pair correlation increases with temperature: 
in the DC regime temperature has a destructive role on 
the correlation. 



VI. DECOHERENT REGIMES 



In the limit of the weakly interacting Bose gas, (7 ^ 
min{T^ , y/r}) or in the high temperature regime (t <C 
max {1,7^}), the phase and density fluctuations are 
large. Therefore one can notice that in the YYT equa- 
tions [sj the limit of either high temperature, = 6 <^ 
1 with finite c, or low coupling c = 5 <^ 1 with T 7^ 0, 
one recovers the thermodynamics of two ideal Bose gases 
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FIG. 5: Polarization of the 2 component Bose gas as a function of the interaction strength 7 for fixed values of ^^2^^ ) 
and for fixed ratios 



up to 0{5'^). In this limit the convolutions of a function 
g with the kernels described in equations ([3| become: 



dk' 



nc/{2VT) 



7r(nc/(2Vr))2 + (fc-fc')^ 



gikVT) 



dk' 



■g{kWT) 



(16) 



T 



cosh(^VT(fc-A:')) 



•g(fcVT) 



cff(fcVT) 



(17) 



Furthermore the Gibbs free energy resulting from this 



simplified system is: 
G T 



L 



+ 



27r 
T 
2^ 
T 
2ti 



In 



1 + e-^(^-)/^ 



In 
In 



1 _ g(Mi-fe')/T 
_ g(A'2-fe')/T 



dk 
dk 



(18) 



which is the sum of the Gibbs energy of two ideal Bose 
gases. 

First order corrections can then be effectively described 
using perturbation theory and the reduced temperature, 
r, allows one to distinguish between the de- 



(711+712) 

coherent quantum regime (DQ) for ^7 ^ r ^ 1 and the 
decoherent classical (DC) regime with r ^ max{l,7'^} 

I1I1I171[5S|. 

We use Feynman diagrams to express the perturbed 
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Gibbs free energy. An explicit expression is then calcu- 
lated for the local pair correlation function, 5^^) in the 
two decoherent regimes (DQ & DC) to first order. 



m 



The partition function of the 2 component ID Bose gas 
the Feynman path integral formalism is 



FIG. 6: The local pair correlation g'^' of the 2 component 
Bose gas as a function of the interaction strength 7 and for 
fixed values of Q and for fixed ratio 
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FIG. 7; Comparison between the numerical (thick lines) and 
the analytical results (thin lines) in the decoherent quantum 
regime, ^ r <C 1 (|2 2|| a nd in the decoherent classical 
regime, r ^ max{l,7^} ( |23[ ). The reduced temperature is 
fixed either to r = 0.1 with — 0.1, 1 for a decoherent quan- 
tum gas or to r = 1000 with Q, — 1, 10'' for a classical de- 
coherent gas. At high value of the interaction strength when 
7 2> 1, the 2CBG enter a fermionization regime. 



Z = y I?(*,*)e-^[*'*l 
S-i*,*] = j^'^^^y dr J2 ^ adr'f a -ni'i?,'^) (19) 



where 'S'(r, r) is a space and imaginary time-dependent 
spin- 1/2 field and H is the Hamiltonian density from ([I]). 
At first order, the correction to the Gibbs free energy fol- 
lowing from Wick's theorem is G*^^-' — 2c 



0{c^) with the free linear density of the a-th component: 

n2 — y^i- 1,2^-2 /n — I — and the total free linear 

density: nP — n°. For the second order in c, the dia- 
grammatic representation gives five contributions shown 
in figure |8] that give the free energy density corrections: 



G(2) = -'- 



E 

a=0,l 



V 



0{c^) (20) 



d) 



10 




FIG. 8: Connected Feynman diagrams in second order perturbation of interaction. Labels are related to terms of eq. (20 1 



where the first three terms correspond to the diagrams containing the double polarization bubbles which are de- 
6), a) and d) and where c) and e) provide the last terms fined as 



/CO / pOG /• OO 

dk dlGa.,n+n{k + l)GaAl)Y. 



.,m-\-n' 



(21) 



with the Green function Ga,n{l) 



The 



local pair correlation results from equation ( 10 ) and an 
analytic expression as a function of c, T, /i^ is given in 
the two decoherent regimes. In the DQ regime, ^ <C 
T ^ 1 and /ii , /i2 ^ T, by taking the leading order 
in the Bose occupation number, the free linear density 



is = 2y^^^ ^"^^ double polarization bubble is 

^,6 = ° ° ° '' . For a compact nota- 

tion we define the atli component reduced temperature 
by Ta — -57. We find so the local pair correlation to be 



P, 



■,(2) 



47 



TIT2 



.02 ^r? ^ 



47- 



4( 



1 1 



T1T2 



-0(7') (22) 



In the DC regime, r ^ max{l,7^} and /ii,/i2 S> T, pair correlator becomes 
the bosonic occupation number becomes the Boltzmann 

distribution and = ./s^, Pa.f, = The 



i(2) - 



.02 



47 



.02 



-0(7') (23) 



In order to illustrate this result, we compa re t he v alue 
of 5^^^ computed to the first order in 7 (eqj22] and [23| 
with the numerical results using ( 10 ) in figure [7} The 
curves calculated at fixed particle density and reduced 
temperature show that the numerical results follow nicely 



the analytical expansion until either the thermal fluctu- 
ations become too strong for the DQ gas (7 ^ 10~^) 
or the charge fluctuations become important in the DC 
regime when 7 ~ 10. As the interaction strength in- 
creases, we progressively switch to a high temperature 
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Tonks- Girardeau like fermionization regime for the DC 
curves and to a ferromagnetic fermionization for the DQ 
case. It would be interesting as well to compare the re- 
sults for a DQ gas in very low fl such that the polariza- 
tion is low and 5*^^^ reaches the value 3/2 but this implies 
a calculation for a very high number of functions, Umax 
with a large number of points. We couldn't afford then 
the number of iterations necessary to have a converged 
solution. 



VII. TONKS-GIRARDEAU REGIME 

In the extreme case of impenetrable particles (7 — )• 00), 
M. Girardeau [21] showed the correspondence between 
impenetrable Bose and Fermi wave functions. While the 
statistics of the bosons wave function remains symmetric, 
there is no more overlap between the neighbor particles. 
In the case of 2CBG, the charge part of the wave function 
behave like a one-component free fermion gas and non- 
interacting distinguishable spin— 1/2 since any spin-spin 
exchange vanishes [Ml EH HI- In both ICBG and 2CBG, 
the local density-density correlation function then natu- 
rally vanishes since there is no double space occupancy. 

In the strong coupling regime (7 ^ max{l, v^^)) with 
quantum degeneracy (r ^ 1), the finite-temperature cor- 
rections are markedly different in a Lieb-Liniger gas [46] 
and the spinor Bose gas [34j with a different exponent. 
In figure [9j bottom part, we represent this analytical re- 
sult at zero temperature (thin line) next to numerical 
results with decreasing temperature for fixed density of 
particles. We observe that for 7 ^ 10, the value of g^^^^ 
decreases with t and converges to this T — analyti- 
cal result where the 2CBG is ferromagnetic and doesn't 
depend on f2. For a high-temperature fermionization 
(7^ ^ T ^ 1), Kheruntsyan et al. [IS] give the first 
order correction in r/7^ for g^^^ for a Lieb-Liniger gas. 
However the approach of free fermions with a I/7 per- 
turbation is not applicable in the two component case, 
therefore the correction to the fermionization regime are 
unknown. In figure [O] top part, we show next to the 
1 component asymptotic curve, the decay of g^'^^ for a 
fixed reduced temperature and various relative chemical 
potential. As reaches 100, the 2CBG polarization is 
saturated and the correlator decays like a Lieb-Liniger 
gas. 



VIII. CONCLUSION 

In conclusion, we have studied the equilibrium thermo- 
dynamic properties of exactly solvable interacting one- 
dimensional two-component Bose gas systems as a func- 
tion of their external canonical or grand canonical param- 
eters (either temperature, interaction strength and total 




10 y 100 



FIG. 9: In the Tonks-Girardeau regime, the pair-pair corre- 
lation decay to zero when the correlation strength becomes 
large (7 ^ max(l, \/t))- The top graph shows the high- 
temperature fermionization and its dependence in Q. The 
ICBG asymptotic behavior in r/7^ is also shown in compar- 
ison with the fully polarized 2CBG {Q ^ T). In the bottom 
graph we show curves with different reduced temperatures 
that we compare with the analytical expression at zero tem- 
perature [H3]. The 2CBG is here ferromagnetic and the value 
of g{2) doesn't depend on the relative chemical potential. 



and relative chemical potential or temperature, interac- 
tion strength, densiy of particle and relative chemical po- 
tential). Our method was based on the solution of ther- 
modynamic Bethe ansatz equations and yields quantita- 
tive predictions which should be experimentally accessi- 
ble using cold atomic systems. We particularly would 
like to clarify that solving the non linear integrable equa- 
tions is possible with a very good control of numerical 
precision. 

Note: as our manuscript was being completed, a dif- 
ferent but equivalent set of equations was proposed in 
[S3]. While this set of equations is at first sight more 
economical, we find and demonstrate here that the solu- 
tion of the infinite set of TBA equations is feasible and 
practical, robust and reliable. The TBA dressed energies 
in [3] are relatively smooth functions of a real variable, 
while the functions of .54J are are of a complex variable. 
The computational effect required by the two methods 
are thus probably comparable. On the other hand, the 
fact that results from this alternate method coincide with 
our results here (and our earlier summary [lU]) interest- 
ingly confirms that the string hypothesis can be trusted 
when computing equilibrium thermodynamic results, as 
expected from general arguments based on the structure 
of the Bethe equations [55) . 
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Appendix A: Thermodynamics from Bethe Ansatz 

We model the system of 2 component bosons with 
SU(2) bosonic fields evolving in a 1 dimensional contin- 
uum space of length L with a delta-function interaction. 
The Hamiltonian is then: 



^° ae{-l,l} 
6,aG{-ia} 



(Al) 



With c 



^'^^2™ , 5iD the ID coupling constant and m 
the mass of the bosons. 



A" 



Integrating the string structure, j 
1 — 2j), j — 1, . . . , n, in the scattering equations in ^ 
and defining e„(A) = A+len/2 ' scattering equations 
become: 



fcp-Ag+'f ^ „ -pr r ei(A)el(A)eL„.+4(A) • . • eL-2(A)e2„(A), m = n 

-^^fcp-Ag-i^ ^ ' i-Ll e„_,^(A)e2_ (A)e2 4(A)...e2 2(A)e„+„(A),m7^n 

p— 1 K ^ m,p 

(A2) 

with the notation: A = AJJ - A™. 

In logarithm form, with i ln(e„(A)) = (/!'ti(A) = — tt + 2 atan( j^^A^), we have: 



^3 

N 



1=1 n=la=l 



p=l 



1 ^ ^ J 2</.2(A2 - A-) + 204(A2 - A^) . . . (/.2„(AS - A^), m = n 

^ „^i^i I -^In-mKAS - A^) + 20|„_„|+2(AS - A^') . . . (AS - A^), m ^ n 

(A3) 



r 



{Ij} is a set of N numbers in K + ^ and {J^} are M„ sets 
of n numbers in K (K -|- i) if M„ is even (odd). These 
Bethe equations map the sets {Ij},{Ja} to the set of 
rapidities and isospin-rapidities, {kj} and{AJJ}. 



1. Thermodynamic Hmit 



L ^ 

j 

1 in 



1 7" 

^"(A') = 7E'^(A-A^(^),Vn . (A4) 



In the limit N, L oo with the ratio ^ kept con- 
stant, the sets of rapidities {{kj} and {A^}) and quan- 
tum numbers {{Ij} and {J^}) are replaced by continu- 
ous functions of particle root densities in real parameter 
space: 



Hole densities, ph-, are similarly defined from the com- 
plementary sets {Ii},{Ja}: and the total root densities 
are pt{k) = ph{k) + ph{k) and a," (A) = al{K) + a" (A). 
The thermodynamic limit allows one to replace the dis- 
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Crete sum by an integral over a continuum. 



Pt{x) = ^ 



^ L Th.L. 



namic limit is the assumption that all the density func- 
tions are in C°° . However some sets of {/j}, {>/"} that 



dx'5{x 



Th.L. 



(A5) 



and the indexations of the rapidities by the quantum 
numbers, {kj{j^),K'^{-f^)) become continuous functions 
: k{x) and A"'(y). An important point of the thcrmody- 



are solutions of ( A2 ), could provide no differentiable func- 



tions. For instance if all the rapidities are grouped in a 
block Fermi-sea like. But the role of these solutions play 
a negligible role in the thermodynamic limit due to the 
fact that their weight in the set of all solutions goes to 
zero. Physically, they represent the solutions with low 
entropy. Finally, under the thermodynamic limit ( A3 ) 
becomes: 



k{x) 

(pnik - A)p{k)dk 



/CO ^ /"OO 

Mk - k')p{k')dk' + E / 't^mik - A)a"'{A)dA 

2Try 

r ,w/ 202(A-A') + ...02„(A-A'), m = n 

hj~oo W|„-h(A-A') + ...0„+„.(A-A'), m^n 



2. YYT equations 

Following the method of C. N. Yang and C. P. Yang 
3] , the equilibrium state is determined by minimization 



(A6) 



r 



of the Gibbs free energy in the grand canonical ensemble. 
With G the Gibbs free energy, E the internal energy, S 
the entropy, we have: 



G = E-TS - piNi- H2N2 

- = / dkp{k)e 

S 
L 



Hini + p2n2 



/oo ^ poo 

dk [{p + p^) Hp + py,) - pHp) - p^ in(p;,)] + Y. [(^" + O + O - - 

j dkn{^-2'Y^ ncj" j -f /ip (A7) 



with L the length of our system, rii = ^ the density of condition of equilibrium is then: 
ith component particles and p = ta^Hz. , O = Bi^Hi. , The 



, dG BG 

Op-^ h Op/i- 

dp dph 



dG ^ ^dG 



(A8) 



p.cr^solutioii of BE 
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from which one derives: 

e(fc) = fc2-/i-17-r-a2* ln(l + e^'l^) - ^ T • a„ * ln(l + e"^"/^) 
en{k) = 2nn + T- a„ * ln(l + e"^/'^) + T • ^ * ln(l + e-^"/^), n = 1, 2, . . . (A9) 



with 





f 2a2(A) + 2a4(A) . 


..a2„(A), TO = n 


\ a\n~m\{^) + 2a|„_ 


-m|+2(A) . . .a„+,„(A),TO ^ n 


e(fc) 


/5(fc) 




e„(fc) 







(AlO) 
(All) 
(A12) 



The term Tm„ which implies a couphng between every the system is partiaUy decoupled and this term disap- 
Sm{k) would severely slow down any numerical solving. pears: 
But following the development of M. Takahashi 144j . 



£(fe) 
£„(fc) 
£l(fc) 



* In 



e-H-n-T-(^a2*ln [l + exp(-|;)]) (fc)-r^ (a„ 

71=1 

{/ * (in [l + exp(^)] + In [l + exp(^)] ) } (k), {n ^ 1) 
{/ * (in [l + exp(^)] + In [l + exp(-|;)] ) } (fc) 



1 + exp( 



T 



(fc) 



T 
2c 



T 



(A13) 
(A14) 
(A15) 



with the convolution notation: (/ * = / /(fc — 

k')g{k')dk', a„{k) = f{k) = l/cosh(f fc). 

We can easly calculate the two asymptotic limit 

similarly to the results for the istropic spin chain of 



M.Takahashi 



lim ^ 

n—^oo Ji 



hm e„(fc)(^ £-) 



= 2n 



2VLn + T • In 



l-exp(-^(n + l)) 



1 — cxp(- 



' T ' 



exp(- 



2VL 

) 

T 



(A16) 
(A17) 



Appendix B: Dressed energy derivatives 



The derivatives of e{k) and £„(fc) are useful for calculation of free energy derivatives. Differentiating ([s]) and ( A17) 
hy 1/ = fi, we get: 
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dv 



(k) = 



lim 

n— >oo 



dv 

dsi 
dv 

den/dn{k) 
n 



-ik) 
-(fc) 



-1 + ^02 * 

/* 

/* 



1 

2c 
1 

2c 



de/dv 
l + cxp(f; 



^ n=l ^ 



den/dv 
+ exp(f-) 



+ 



d£n-i/dv 



l + exp(-£^) l + exp(-£^) 
/dv de I d\x 



1 + exp(- 



■^1 



1 + exp(|;) 



(fc) 



2(l-exp(-f (n+1))) 



(1 - exp(-^)) (1 - exp(-^n)) (l - exp(- 



f(n + 2))) 



2f2 2f2 2f7 

n - n • exp(- — (n + 2)) + (n + 2) • exp(- — (n + 1)) - (n + 2) • exp(- — ) 



(Bl) 



with the derivatives of the asymptotes for = being identically zero. Differentiating by T gives: 



5£ 

den 
dT 

dsi 



(k) 
(k) 
(fc) 



lim 

n— >-oo 



dT 
den/dT{k) 



dT 



e{k) -k'^ + fi + fl 



02 



T 2c 

gl(fc) , 1 

T 2c 



, , den+i/dT - En+i/T den-i/dT - e^-i/T 

J * I — / e„^iN H 



l + exp(-^) ' l + exp( 

de2/dT - e-ilT _ de/dT-e/T 
1 +exp(-f ) l + exp(f) 



En-l ' 

T - 



exp(^) 

(fc), (n 9^ 1) 



(fc) 



= 



In 



-exp(-f^(n + l)) 
l_exp(-f ) 



+ 



217 ■ 



2 e-P(- f ) ^"'IZttt))^ - „ exp(- f n) - 2(n + 1) exp(- f (n + 1)) 



exp(-^(n+l)) 



(l_exp(-3|l))^ 



( 



l-exp(-gp («+!)) \ _ 
l-exp(-^^) I ^ 



(B2) 
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And by c: 



lim 



deik) 
dc 



denjk) 
dc 



dei{k) 
dc 



den/dc{k) 



n—^oo ji 



dc{k) 



oo 

l + cxp(-|,)])(A:)-T^ 



'dan 
dc 



^ de/d, 
y l+exp(f) 



^ n=l ^ 

In 



1 + exp(-^; 



T 


'(df 


Yc 


\dc 


1 


/* ( 




Yc 


T 


'(df 


Yc 


\dc 


1 






Yc 



I'" 

d£n+l/dc 



1 +cxp(^^) 

1 + exp(-2^) 



-(^) 



1 + cxp( — 
(fc),(n^l 



{k) 



In 



l+cxp(|) 



-In 



/ de2/dc dejdc \ 

lvl + exp(-f ) ~ l + exp(f.)j 



1 + exp(-- 
(fc) 



(k) 








(fc) 



(B3) 



Appendix C: Covariance under the parameter c This allows one to reduce our parameter space to 

{/z, ri, r} and put c = 1 by default. 

The thermodynamics of the system depend on four pa- 
rameters: (c, ^, ri,T). Or we can easly deduce form ([S]) 
that they are covariant under renormalisation by c, i.e.: 

c c c 

ni{c,n,fl,T) = crii(l, ^, ^, = 1,2 (CI) 
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